Examine fluxes

# Define functions to convert to pM or to #/cell
number_per_cell <- function(conc, tissue){
  ESAV = ifelse(tissue == "main", 73, 105)
  conc*6.0230e+23*1e-5/ESAV;
}

convert_to_pmol <- function(conc, tissue){
  K_AV = ifelse(tissue == "primary", 0.58, 0.03)
  conc/K_AV*1e15
}

convert_secretion_mol_per_cm3_tissue <- function(q_V165, tissue = "primary"){
  q_V165  * 1534 / (1e-5 * 6.023e23)
}

quick_bar <- function(df, m,convert = "#/cell", tissue = "primary", position = "stack"){
  converter = switch(convert, "#/cell" = number_per_cell, "pM" = convert_to_pmol)
  y_lab = ifelse(position == "fill", "proportion", convert)
  
  df <- df %>% 
    filter(compartment == tissue, grepl(m, molecule), q_V165 < 1e4) %>% 
    group_by(q_V165) %>% 
    filter(time == max(time)) %>% 
    ungroup
  
  if(m == "V165"){
    df <- df %>% pivot_wider(names_from = molecule, values_from = concentration) %>%  
      group_by(q_V165) %>% 
      summarise(Free = V165, 
                `Bound to M` = sum(across(starts_with("M") & !contains("R"))), 
                MVR = sum(across(starts_with("Mebm_V165_"))),
                RVN = sum(across(starts_with("R2_V165_"))),
                `Receptor Bound` = R1_V165 + R2_V165 + N1_V165 + N2_V165) %>% 
      pivot_longer(!q_V165, names_to = "molecule", values_to = "concentration")
  }
  
  df %>% 
    mutate(value = converter(concentration, tissue)) %>% 
    ggplot(aes("", value, fill = molecule)) + 
    geom_bar(stat = "identity", position = position) + 
    labs(x = "", y = y_lab, fill = "Complex",
         title = str_interp("${m} in ${tissue}")) + 
    facet_wrap(~q_V165, scales = "free_y") +
    theme(legend.position = c(1,0), legend.justification = c(1,0))
}

Fluxes into and out of the primary tumor compartment.

Using the original model parameters, I calculated secretion, lymphatic drainage, internalization, and permeability at steady state. Fluxes were normalized by dividing by the rate of secretion (ie. a value of 1 indicates that it is equal to the rate of secretion). I found that with the original model parameters, lymphatic drainage was almost equal to secretion unless the rate of R2 production was very large. This indicates that the majority of VEGF being secreted is carried out of the tissue without the chance to bind, which is not what we want.

fluxes <- read_csv(here("debug_results", "calculate_fluxes.csv")) %>%
  mutate(q_V165 = factor(q_V165), kprod_R2 = factor(kprod_R2))

# Plotting function
plot_norm_flux <- function(flux){
  flux %>%
    pivot_longer(4:6) %>%
    mutate(norm_value = abs(value/Secretion)) %>%
    ggplot(aes(kprod_R2, norm_value, color = name)) +
    geom_point(size = 3) +
    facet_wrap(~q_V165, labeller= "label_both")  +
    labs(x = "Production of R2 (moles/cm^3 tissue)", y = "Normalized Flux", color = "Flux") +
    theme(legend.position = c(1,0), legend.justification = c(1,0))
}


plot_norm_flux(fluxes)

I updated the lymphatic drainage for the primary tumor as \(5.6\cdot10^{-7}\cdot 5\cdot 6.4=1.79\cdot10^{-5} cm^3/s\) according to lymphatic drainage when awake multiplied by the tissue mass. Now, internalization dominates unless V165 secretion is greater than the physiological limits, which is what we expect.

new_param_fluxes <- read_csv(here("debug_results", "calculate_fluxes_only_R2_new_lymph.csv")) %>%
  mutate(q_V165 = factor(q_V165), kprod_R2 = factor(kprod_R2))

plot_norm_flux(new_param_fluxes) +
  ggtitle("Fluxes with new model parameters")

orig_simulations <- read_csv(here("debug_results", "simulation_results_debug_R2.csv"))
new_simulations <- read_csv(here("debug_results", "simulation_results_debug_R2_only_R2_new_lymph.csv"))

orig_simulations_long <- orig_simulations %>%
  pivot_longer(!c(q_V165, kprod_R2, time), values_to = "concentration",
               names_to = c("compartment", "molecule"), names_sep = "[.]")
new_simulations_long <- new_simulations %>%
  pivot_longer(!c(q_V165, kprod_R2, time), values_to = "concentration",
               names_to = c("compartment", "molecule"), names_sep = "[.]")

Concentration profile of V165

The concentration of free V165 with the baseline secretion and R2 production increased from 0.3 to >4 pM. In addition, there is a greater decrease in V165 concentration as R2 production increases.

orig_free_V165 <- orig_simulations %>%
  ggplot(aes(time, convert_to_pmol(primary.V165, "primary"), color = as.factor(kprod_R2))) +
  geom_line() +
  facet_wrap(~q_V165, labeller = "label_both", scales = "free_y", nrow = 1) +
  ggtitle("Original Parameters")

new_free_V165 <- new_simulations %>%
  ggplot(aes(time, convert_to_pmol(primary.V165, "primary"), color = as.factor(kprod_R2))) +
  geom_line() +
  facet_wrap(~q_V165, labeller = "label_both", scales = "free_y", nrow = 1) +
  ggtitle("Updated Parameters")

orig_free_V165 + new_free_V165 +
  plot_layout(nrow = 2, guides= "collect") &
  labs(y = "Free V165 (pM)", x = "time (s)", color = "R2 Production Rate")

quick_bar <- function(df_long, m, tissue = "primary", pos = "stack", convert = "pM", x = NULL, scales = "free_y", outline = NULL){
  if(is.null(x)) x <- ifelse(m == "V165", "kprod_R2", "q_V165")
  facet = c("kprod_R2", "q_V165")[c("kprod_R2", "q_V165") != x]
  
  x_lab = switch(x, "kprod_R2" = "R2 Production (moles/cm^3 tissue)", "q_V165" = "V165 Secretion (#/cell)")
  y_lab = switch(pos, "stack" = str_interp("${m} (${convert})"), "fill" = str_interp("Fraction of ${m}"))
  converter = switch(convert, "pM" = convert_to_pmol, "#/cell" = number_per_cell)
  
  p <- df_long %>%
    group_by(q_V165, kprod_R2) %>%
    filter(time == max(time), grepl(m, molecule), compartment == tissue) %>%
    ggplot(aes_string(str_interp("as.factor(${x})"), "converter(concentration, tissue)", fill = "molecule")) +
    facet_wrap(as.formula(paste0("~", facet)), scales = scales, nrow = 1, labeller = "label_both") +
    labs(x = x_lab, y = y_lab)
  
  if(!is.null(outline)){
    p + geom_bar(stat = "identity", pos = pos, aes(color = grepl(outline, molecule)))  +
      scale_color_manual(values = c("TRUE" = "black", "FALSE" = NULL)) +
      labs(color = "R2 Complex")
  }else{
    p + geom_bar(stat = "identity", pos = pos)
  }
}

The breakdown of total V165 concentration is as follows, with the black outline denoting V165 in a complex with R2. There is more V164 total now.

quick_bar(orig_simulations_long, "V165", outline = "R2") +
  labs(x = "", title = "Original Parameters") +
  quick_bar(new_simulations_long, "V165", outline = "R2") +
  labs( title = "Updated Parameters") +
  plot_layout(nrow = 2, guides = "collect") &
  theme(panel.grid.major.y = element_line())

quick_bar(orig_simulations_long %>% filter(molecule != "V165"), "V165", outline = "R2") +
  labs(x = "", title = "Original Parameters") +
  quick_bar(new_simulations_long  %>% filter(molecule != "V165"), "V165", outline = "R2") +
  labs( title = "Updated Parameters") +
  plot_layout(nrow = 2, guides = "collect") &
  theme(panel.grid.major.y = element_line()) &
  labs(y = "Bound V165 (pM)")

I also looked at the fraction of V165 in each complex. I spent a lot of time thinking about why the fraction of bound V165 remains the same for the lower values of V165 secretion. I found that in these cases, the amount of free and bound V165 increases proportionally to the increase in total V165 so the fraction stays the same. I don’t know if this makes sense biologically because we have more binding and therefore more internalization of V165.

quick_bar(orig_simulations_long, "V165", outline = "R2", pos = "fill") +
  labs(x = "", title = "Original Parameters") +
  quick_bar(new_simulations_long, "V165", outline = "R2", pos = "fill") +
  labs( title = "Updated Parameters") +
  plot_layout(nrow = 2, guides = "collect") &
  theme(panel.grid.major.y = element_line())

Compare simulations

Below is a direct comparison of the concentrations of free, bound, and total V165 between the two simulations. Across all parameter combinations, both the total amount of V165 and free V165 is increased; more so when V165 secretion is greater than R2 production. We see a similar trend in the increase in receptor and matrix-bound V165. The bottom left corner of the plot for receptor-bound V165 is very interesting.

plot_df <- merge(orig_simulations_long, new_simulations_long, by = 1:5) %>% 
  filter(time == max(time), grepl("V165", molecule), compartment == "primary") %>% 
  pivot_longer(starts_with("conc")) %>%
  # mutate(category = ifelse(molecule == "V165", "Free V165", 
  #                          ifelse(molecule %in% paste(c("Mebm", "Mecm", "Mpbm"), "V165", sep = "_"), 
  #                                 "Matrix-bound V165", "Receptor-bound V165")),
  #        value = convert_to_pmol(value, "primary")) %>% 
  mutate(category = ifelse(molecule == "V165", "Free V165", 
                           ifelse(molecule %in% paste(c("Mebm", "Mecm", "Mpbm"), "V165", sep = "_"), 
                                  "Matrix-bound V165", 
                                  ifelse(grepl("Mebm", molecule) & grepl("R2", molecule), "MVR",
                                         ifelse(grepl("R2", molecule), "R2-bound V165",
                                                "Other Receptor-bound V165")))),
         value = convert_to_pmol(value, "primary")) %>% 
  group_by(q_V165, kprod_R2, name, category) %>%
  summarise(value = sum(value)) %>% 
  pivot_wider()  

totals <- plot_df %>% summarise(concentration.x = sum(concentration.x), 
                                concentration.y = sum(concentration.y)) %>% 
  mutate(category = "Total V165")

plot_df %>% 
  ggplot(aes(concentration.x, concentration.y, 
             color = as.factor(q_V165), shape = as.factor(kprod_R2))) + 
  geom_point() + 
  geom_point(data = totals ) +
  geom_abline() + 
  facet_wrap(~category, scales = "free") + 
  scale_x_log10() + scale_y_log10() + 
  labs(x = "Old Paramaeters", y= "New Parameters", 
       color = "V165 secretion\n(#/cell)", shape = "R2 Production\n(moles/cm^3 tissue)",
       title = "Compare concentrations") 

plot_df %>% 
  ggplot(aes(concentration.x, concentration.y, 
             color = log(kprod_R2/convert_secretion_mol_per_cm3_tissue(q_V165)))) + 
  geom_point() + 
  geom_point(data = totals ) +
  geom_abline() + 
  facet_wrap(~category, scales = "free") + 
  scale_color_gradient2() + 
  scale_x_log10() + scale_y_log10() +
  labs(x = "Old Paramaeters", y= "New Parameters", 
       color = "log(kprod_R2/q_V165)",
       title = "Compare concentrations") 

I also compared the fraction of total V165 in each state: free, matrix-bound or receptor-bound. I found that when R2 production is greater than V165 secretion, these values do not differ between the original and new simulations.

plot_df %>% 
  mutate(concentration.x = concentration.x/sum(concentration.x),
         concentration.y = concentration.y/sum(concentration.y)) %>%
  ggplot(aes(concentration.x, concentration.y, 
             color = as.factor(q_V165), shape = as.factor(kprod_R2))) + 
  geom_point() + 
  geom_abline() + 
  facet_wrap(~category) + 
  coord_fixed() +
  labs(x = "Old Paramaeters", y= "New Parameters", 
       color = "V165 secretion\n(#/cell)", shape = "R2 Production\n(moles/cm^3 tissue)",
       title = "Compare concentration fractions") 

plot_df %>% 
  mutate(concentration.x = concentration.x/sum(concentration.x),
         concentration.y = concentration.y/sum(concentration.y)) %>%
  ggplot(aes(concentration.x, concentration.y, 
             color = log(kprod_R2/convert_secretion_mol_per_cm3_tissue(q_V165)))) + 
  geom_point() + 
  geom_abline() + 
  facet_wrap(~category) + 
  labs(color = "log(kprod_R2/q_V165)") +
  scale_color_gradient2() + coord_fixed() + 
  labs(x = "Old Paramaeters", y= "New Parameters", 
       color = "log(kprod_R2/q_V165)",
       title = "Compare concentration fractions") 

Concentration profile of R2

Bound R2

We see an increase in bound R2 compared to the original model parameters, as we would expect given the increased internalization flux.

orig_bound_R2 <- orig_simulations %>%
  rowwise %>%
  mutate(bound_R2 = sum(c_across(starts_with("primary") & contains("R2") & contains("_")))
  ) %>%
  ggplot(aes(time, number_per_cell(bound_R2, "primary"), color = as.factor(kprod_R2))) +
  geom_line() +
  facet_wrap(~q_V165, labeller = "label_both", scales = "free_y", nrow = 1)

new_bound_R2 <- new_simulations %>%
  rowwise %>%
  mutate(bound_R2 = sum(c_across(starts_with("primary") & contains("R2") & contains("_")))
  ) %>%
  ggplot(aes(time, number_per_cell(bound_R2, "primary"), color = as.factor(kprod_R2))) +
  geom_line() +
  facet_wrap(~q_V165, labeller = "label_both", scales = "free_y", nrow = 1)

orig_bound_R2 + new_bound_R2 +
  plot_layout(nrow = 2, guides= "collect") &
  labs(y = "Bound R2 (#/cell)", x = "time (s)", color = "R2 Production Rate")

As secretion increases, the total amount of R2 decreases due to increased binding and internalization as we would expect.

quick_bar(orig_simulations_long, "R2", convert = "#/cell") +
  labs(x = "", title = "Original Parameters") +
  quick_bar(new_simulations_long, "R2", convert = "#/cell") +
  labs( title = "Updated Parameters") +
  plot_layout(nrow = 2, guides = "collect") &
  theme(panel.grid.major.y = element_line())

With the new parameters, there is more bound R2 both in concentration and fraction of total R2 at the lower secretion levels. This is because increased binding leads to increased internalization and less total R2 so the amount of bound R2 becomes a greater fraction.

quick_bar(orig_simulations_long %>% filter(molecule != "R2"), "R2", convert = "#/cell") +
  labs(x = "", title = "Original Parameters") +
  quick_bar(new_simulations_long  %>% filter(molecule != "R2"), "R2", convert = "#/cell") +
  labs( title = "Updated Parameters") +
  plot_layout(nrow = 2, guides = "collect") &
  theme(panel.grid.major.y = element_line()) &
  labs(y = "Bound R2 (#/cell)")

quick_bar(orig_simulations_long, "R2", convert = "#/cell", pos = "fill") +
  labs(x = "", title = "Original Parameters") +
  quick_bar(new_simulations_long, "R2", convert = "#/cell", pos = "fill") +
  labs( title = "Updated Parameters") +
  plot_layout(nrow = 2, guides = "collect") &
  theme(panel.grid.major.y = element_line())

Free R2

With the new model parameters, the amount of free R2 behaves more closely to the model with no permeability/lymphatic drainage when VEGF secretion and R2 production are large.

orig_free_R2 <- orig_simulations %>%
  ggplot(aes(time, number_per_cell(primary.R2, "primary"), color = as.factor(kprod_R2))) +
  geom_line() +
  facet_wrap(~q_V165, labeller = "label_both", scales = "free_y", nrow = 1)

new_free_R2 <- new_simulations %>%
  ggplot(aes(time, number_per_cell(primary.R2, "primary"), color = as.factor(kprod_R2))) +
  geom_line() +
  facet_wrap(~q_V165, labeller = "label_both", scales = "free_y", nrow = 1)

orig_free_R2 + new_free_R2 +
  plot_layout(nrow = 2, guides= "collect") &
  labs(y = "Free R2 (#/cell)", x = "time (s)", color = "R2 Production Rate")

Compare levels of Free Ligand to Bound Receptors

I compared the amount of free ligand to receptor-bound ligand at every timepoint from the model with updated parameters.

# pdf(here("fig", "free_vs_bound_ligand.pdf"), width = 12, height = 12)
compartment <- "primary"
new_simulations %>% 
  select(kprod_R2, q_V165, time, starts_with("primary")) %>% rowwise %>% 
  mutate(bound_R = sum(c_across(starts_with(compartment) & contains("V165") & (contains("R1") | contains("R2") | contains("N1") | contains("N2")))),
         total_R = sum(c_across(starts_with(compartment) & (contains("R1") | contains("R2") | contains("N1") | contains("N2")))),
         total_V165 = sum(c_across(starts_with(compartment) & contains("V165"))),
         free_V165 = c_across(matches(paste0(compartment, ".V165")))) %>% 
  ggplot(aes(convert_to_pmol(free_V165, "primary"), convert_to_pmol(bound_R, "primary"))) + 
  geom_point(aes(color = time)) +
  facet_wrap(kprod_R2 ~ q_V165, scales = "free", labeller = "label_both") +
  labs(x = "Free V165 (pM)", y = "Receptor-bound V165 (pM)")

# dev.off()

I also plotted the steady state values of free and receptor-bound V165 from all the simulations I’ve run. I plotted both the total concentration of these quantities and the fraction of receptor-bound V165 vs the fraction of free V165. Neither of these fractions include matrix-bound V165, so they are not expected to sum to one. I colored them both by V165 secretion and the log of the ratio of R2 production and V165 secretion such that positive numbers represent greater R2 Production than V165 secretion and vice versa.

plot_R_v_free_ligand <- function(df, prop = F, R = "Receptor", compartment = "primary", scale = "free", shape = "kprod_R2",color = "q_V165"){
  # R denotes receptor to plot on x-axis
  
  if(R == "Receptor"){
    # plot all bound receptors
    df <- df %>%
      mutate(kprod_R2 = factor(kprod_R2), q_V165 = factor(q_V165)) %>%
      group_by(kprod_R2, q_V165) %>%
      filter(time == max(time)) %>%
      rowwise %>%
      summarise(bound_R = sum(c_across(starts_with(compartment) & contains("V165") & (contains("R1") | contains("R2") | contains("N1") | contains("N2")))),
                total_R = sum(c_across(starts_with(compartment) & (contains("R1") | contains("R2") | contains("N1") | contains("N2")))),
                total_V165 = sum(c_across(starts_with(compartment) & contains("V165"))),
                free_V165 = c_across(matches(paste0(compartment, ".V165")))
      ) %>%
      mutate(convert_to_pmol(across(where(is.numeric)), compartment))
  }else{
    df <- df %>%
      mutate(kprod_R2 = factor(kprod_R2), q_V165 = factor(q_V165)) %>%
      group_by(kprod_R2, q_V165) %>%
      filter(time == max(time)) %>%
      rowwise %>%
      summarise(bound_R = sum(c_across(starts_with(compartment) & contains(R) & contains("_V165"))),
                total_R = sum(c_across(starts_with(compartment) & contains(R))),
                total_V165 = sum(c_across(starts_with(compartment) & contains("V165"))),
                free_V165 = c_across(matches(paste0(compartment, ".V165")))
      ) %>%
      mutate(convert_to_pmol(across(where(is.numeric)), compartment))
  }
  
  if(prop){
    # p <- df %>% ggplot(aes(bound_R/total_R, free_V165/total_V165)) +
    #   labs(x = str_interp("Fraction of bound ${R}"), y = "Fraction of Free V165") +
    #   facet_wrap(~q_V165,  labeller = "label_both", nrow = 1, scales = scale)
    p <- df %>% ggplot(aes(y = bound_R/total_V165, x = free_V165/total_V165)) +
      labs(y = str_interp("Fraction of ${R}-bound V165"), x = "Fraction of Free V165") #+
    # facet_wrap(~q_V165,  labeller = "label_both", nrow = 1, scales = scale)
  }else{
    p <- df %>% ggplot(aes(y = bound_R, x = free_V165)) +
      labs(y = str_interp("${R}-bound V165 (pM)"), x = "Free V165 (pM)") #+
    # facet_wrap(~q_V165, labeller = "label_both", nrow = 1, scales = scale)
  }
  
  p + geom_point(aes_string(shape = shape, color = color), size = 2)
}

Free VEGF versus VEGF bound to all receptors

plot_R_v_free_ligand(orig_simulations, F) +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line()) &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical")  &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000))  

plot_R_v_free_ligand(orig_simulations, T) +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical")  &
  coord_fixed() &
  xlim(0,0.6) & ylim(0,0.6)

color <- "log(as.numeric(as.character(kprod_R2))/convert_secretion_mol_per_cm3_tissue(as.numeric(as.character(q_V165))))"


plot_R_v_free_ligand(orig_simulations, F,
                     shape = NULL, color = color) +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F,
                       shape = NULL, color = color) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line()) &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical")  &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000)) &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

plot_R_v_free_ligand(orig_simulations, T, color = color, shape = NULL) +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, color = color, shape = NULL) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical")  &
  coord_fixed() &
  xlim(0,0.6) & ylim(0,0.6) &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

Free V165 Versus V165 bound to R2

Within the same secretion value, increased V165-R2 formation results in less free V165. The fraction of R2-bound V165 generally decreases as free V165 increases, but there are a few parameter combinations that result in ~40% of V165 free and ~0% bound to R2.

plot_R_v_free_ligand(orig_simulations, F, R= "R2") +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "R2") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000)) 

plot_R_v_free_ligand(orig_simulations, T, R= "R2") +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "R2") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  coord_fixed() &
  xlim(0,0.6) & ylim(0,0.6)

plot_R_v_free_ligand(orig_simulations, F, R= "R2",
                     color = color, shape = NULL) +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "R2",
                       color = color, shape = NULL) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000)) &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

plot_R_v_free_ligand(orig_simulations, T, R= "R2",
                     color = color, shape = NULL) +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "R2",
                       color = color, shape = NULL) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  coord_fixed() &
  xlim(0,0.6) & ylim(0,0.6) &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

Free V165 Versus V165 bound to R1

As free V165 increases, so does the amount of R1-bound V165. A very small proportion of V165 is bound to R1 at all simulation values. This is because R1 is mostly bound to neuropilin and V165 does not bind to R1-N.

plot_R_v_free_ligand(orig_simulations, F, R= "R1") +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "R1") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000)) 

plot_R_v_free_ligand(orig_simulations, T, R= "R1") +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "R1") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "veritcal") &
  xlim(0,0.6) & ylim(0,0.004)

plot_R_v_free_ligand(orig_simulations, F, R= "R1",
                     color = color, shape = NULL) +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "R1",
                       color = color, shape = NULL) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000)) &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

plot_R_v_free_ligand(orig_simulations, T, R= "R1",
                     color = color, shape = NULL) +
  labs(title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "R1",
                       color = color, shape = NULL) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "veritcal") &
  xlim(0,0.6) & ylim(0,0.004) &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

Free V165 Versus V165 bound to N1

Binding of V165 to N1 increases as free V165 increases. The fraction of neuropilin-bound N1 decreases as free V165 increases, as we would expect.

plot_R_v_free_ligand(orig_simulations, F, R= "N1") +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "N1") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000)) 

plot_R_v_free_ligand(orig_simulations, T, R= "N1") +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "N1") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  xlim(0,0.6) #& ylim(0,0.004)

plot_R_v_free_ligand(orig_simulations, F, R= "N1",
                     shape = NULL, color = color) +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "N1",
                       shape = NULL, color = color) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000))  &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

plot_R_v_free_ligand(orig_simulations, T, R= "N1",
                     shape = NULL, color = color) +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "N1",
                       shape = NULL, color = color) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  xlim(0,0.6)  &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

Free V165 Versus V165 bound to N2

Binding of V165 to N2 increases as free V165 increases. The fraction of neuropilin-bound N2 decreases as free V165 increases, as we would expect.

plot_R_v_free_ligand(orig_simulations, F, R= "N2") +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "N2") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000)) 

plot_R_v_free_ligand(orig_simulations, T, R= "N2") +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "N2") +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  xlim(0,0.6) #& ylim(0,0.004)

plot_R_v_free_ligand(orig_simulations, F, R= "N2",
                     shape = NULL, color = color) +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, F, R= "N2",
                       shape = NULL, color = color) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  scale_x_log10(breaks = c(0.1, 1, 10, 100, 1000)) & 
  scale_y_log10(breaks = c(0.1, 1, 10, 100, 1000))  &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

plot_R_v_free_ligand(orig_simulations, T, R= "N2",
                     shape = NULL, color = color) +
  labs( title = "Original Parameters") +
  plot_R_v_free_ligand(new_simulations, T, R= "N2",
                       shape = NULL, color = color) +
  ggtitle("Updated Parameters") +
  plot_layout(nrow = 1, guides = "collect") &
  theme(legend.position = "bottom", panel.grid.major.y = element_line(),
        legend.box = "vertical") &
  xlim(0,0.6)  &
  labs(color = "log(kprod_R2/q_V165)") &
  scale_color_gradient2() 

Contour plots

To find the combination of V165 secretion rate and R2 production rate that results in a reasonable amount of total R2 (~3000 receptors/cell) and free V165 (I’m not sure what value we are aiming for here) in the primary compartment, I made contour plots of our updated simulations. The contours for R2 are in terms of #/cell and that of V165 are in terms of pM (Would it be better to plot in terms of pmol/cm^3 tissue?).

contour_df <- new_simulations_long %>% 
  filter(compartment == "primary") %>% 
  group_by(q_V165, kprod_R2) %>% 
  filter(time == max(time)) %>% 
  summarise(total_R2 = number_per_cell(sum(concentration*grepl("R2", molecule)), "primary"),
            free_V165 = convert_to_pmol(sum(concentration*(molecule == "V165")), "primary"),
            perc_bound_R2 = number_per_cell(sum(concentration*(grepl("R2", molecule) & grepl("V165", molecule))), "primary")/total_R2*100)

When using R2 and V165 concentrations directly taken from the model output, the larger values blow up the contours and we do not get any information for the lower quantities. Below, I plotted the points at which we actually ran simulations for reference:

contour_df %>% 
  ungroup %>% 
  # filter(q_V165 !=max(q_V165), kprod_R2 != max(kprod_R2)) %>% 
  ggplot(aes(kprod_R2, q_V165  * 1534 / (1e-5 * 6.023e23))) + 
  geom_contour(aes(z = total_R2, color = after_stat(level)), bins = 25) + 
  scale_color_gradient2("Total R2 (#/cell)", #limits = c(0, 4), 
                        low = "#762A83", mid = "white", high = "#1B7837") +
  # breaks = log(c(3, 30, 3e2, 3e3, 3e4, 3e5)),
  # labels = function(breaks) round(exp(breaks))) +
  
  new_scale("color") +
  geom_contour(aes(z = free_V165, color = after_stat(level)), bins = 25) +
  scale_color_gradient2("Free V165 (pM)",#limits = c(0, 3),
                        low = "#e8cced", high = "#762A83") + #,
  # breaks = log(c(0,1,10,100,1000,10000,100000)),
  # labels = function(breaks) round(exp(breaks))) +
  
  geom_point() +
  labs(x = "Production of R2 (moles/cm^3 tissue/s)", y = "V165 Secretion (moles/cm^3 tissue/s)",
       title = "Contour plot of Total R2 and Free V165")

So, I calculated the contours for V165 and R2 based on the log of their concentration and log-transformed the x and y-axes. The color bars show the concetration values converted back to the original scale (either #/cell or pM).

contour_df %>% 
  ggplot(aes(kprod_R2, q_V165  * 1534 / (1e-5 * 6.023e23))) + 
  geom_point() +
  geom_contour(aes(z = log(total_R2), color = after_stat(level),
                   # Mark 3000 R2:
                   lty = after_stat(level) == after_stat(level)[which.min(abs(after_stat(level) - log(3000)))]),
               bins = 10 ) +
  scale_color_gradient2("Total R2 (#/cell)", 
                        low = "#762A83", mid = "white", high = "#1B7837",#) +
                        breaks = log(c(3, 30, 3e2, 3e3, 3e4, 3e5)),
                        labels = function(breaks) round(exp(breaks))) +
  
  new_scale("color") +
  geom_contour(aes(z = log(free_V165), color = after_stat(level)),
               bins = 10) +
  scale_color_gradient("Free V165 (pM)",
                       low = "#e8cced", high = "#762A83",
                       breaks = log(c(0,1,10,100,1000,10000,100000)),
                       labels = function(breaks) round(exp(breaks))) +
  # geom_point() +
  labs(x = "Production of R2 in moles/cm^3 tissue/s", y = "V165 Secretion in moles/cm^3 tissue/s",
       title = "Contour plot of Total R2 and Free V165 on the log scale", lty = "Total R2 near 3000") +
  scale_x_log10() + scale_y_log10()

The contours intersect and thus we should be able to identify values of q_V165 and kprod_R2 that generate the desired concentrations.

I also looked at the countours for the percent of R2 bound by ligand:

contour_df %>% 
  ggplot(aes(kprod_R2, q_V165  * 1534 / (1e-5 * 6.023e23))) + 
  geom_contour(aes(z = perc_bound_R2, color = after_stat(level)),
               bins = 20) +
  scale_color_gradient2("R2 bound to Ligand (%)",
                        low = "#659473", high = "#1B7837") +#) +
  # breaks = seq(0,100, by = 25)) +
  
  new_scale("color") +
  geom_contour(aes(z = log(free_V165), color = after_stat(level)),
               bins = 10) +
  scale_color_gradient("Free V165 (pM)",
                       low = "#e8cced", high = "#762A83",
                       breaks = log(c(0,1,10,100,1000,10000,100000)),
                       labels = function(breaks) round(exp(breaks))) +
  # geom_point() +
  labs(x = "Production of R2 in moles/cm^3 tissue/s", y = "V165 Secretion in moles/cm^3 tissue/s",
       title = "Contour plot of Total R2 and Free V165 on the log scale") +
  scale_x_log10() + scale_y_log10()

It’s difficult to tell which contours correspond to these desired V165 and R2 concentrations, so I also plotted them individually (I tried to get number labels on the above plot but I was unsuccessful):

contour_df %>% 
  ggplot(aes(kprod_R2, q_V165  * 1534 / (1e-5 * 6.023e23))) + 
  stat_contour(aes(z = log(total_R2), color = as.factor(round(exp(..level..))))) +
  labs(x = "Production of R2 in moles/cm^3 tissue/s", y = "V165 Secretion in moles/cm^3 tissue/s",
       title = "Contour plot of Total R2 and Free V165 on the log scale", lty = "Total R2 near 3000") +
  scale_x_log10() + scale_y_log10() + 
  labs(color = "Total R2 (#/cell)")

contour_df %>% 
  ggplot(aes(kprod_R2, q_V165  * 1534 / (1e-5 * 6.023e23))) + 
  stat_contour(aes(z = log(free_V165), color = as.factor(round(exp(..level..))))) +
  labs(x = "Production of R2 in moles/cm^3 tissue/s", y = "V165 Secretion in moles/cm^3 tissue/s",
       title = "Contour plot of Total R2 and Free V165 on the log scale", 
       lty = "Total R2 near 3000") +
  scale_x_log10() + scale_y_log10() + 
  labs(color = "Free V165 (pM)")